
clear
clear all
clear mata
clear matrix
set more off
set matsize 11000
set maxvar 30000
cap log off
capture log close
set emptycells drop
pause on



***********************************************
* USER DEFINE FILEPATH FOR REPLICATION FOLDER *
***********************************************

global path "C:\Users\wb520443\Dropbox\Research Projects\Global pollution\2_analysis\Replication"



***********************************************


global opts		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))  ) star(* .10 ** .05 *** .01) noobs nocons

global opts_fs		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))) star(* .10 ** .05 *** .01) noobs nocons


global input "$path\data"
global figdat "$path\figdat"
global figures "$path\figures"
global tables "$path\tables"

cap mat define countryfigure=J(105,10,.)
		local row=1
		
	use "$input/rwi_pollution_analysis_collapsed.dta", replace
	
	replace error=1/error
	g both=population*error
	g both2=country_population*error
	
	
	cap encode country, g(country_code)
		label save using "$input/labels.do", replace

	//Mean, mean_ss, and max refer to how PM25 across years and months is collapse (e.g. max is the max pollution level for a RWI point across all years and months in the sample).
	forvalues i=1(1)103{
		
		cap reghdfe mean_pm25 rwi [pweight=both] if country_code==`i', noabsorb
			cap mat def countryfigure[`row',1]=_b[rwi]
			cap mat def countryfigure[`row',2]=_se[rwi]
		cap reghdfe mean_pm25_ss rwi  [pweight=both]if country_code==`i', noabsorb
			cap mat def countryfigure[`row',3]=_b[rwi]
			cap mat def countryfigure[`row',4]=_se[rwi]
		
		local row=`row'+1
	}	
	local row=1
	forvalues i=1(1)103{
		
		cap reghdfe median_pm25 rwi [pweight=both] if country_code==`i', noabsorb
			cap mat def countryfigure_median[`row',1]=_b[rwi]
			cap mat def countryfigure_median[`row',2]=_se[rwi]
		cap reghdfe median_pm25_ss rwi  [pweight=both]if country_code==`i', noabsorb
			cap mat def countryfigure_median[`row',3]=_b[rwi]
			cap mat def countryfigure_median[`row',4]=_se[rwi]
		
		local row=`row'+1
	}	
	local row=1
	forvalues i=1(1)103{
		
		cap reghdfe max_pm25 rwi [pweight=both] if country_code==`i', noabsorb
			cap mat def countryfigure_max[`row',1]=_b[rwi]
			cap mat def countryfigure_max[`row',2]=_se[rwi]
		cap reghdfe max_pm25_ss rwi  [pweight=both]if country_code==`i', noabsorb
			cap mat def countryfigure_max[`row',3]=_b[rwi]
			cap mat def countryfigure_max[`row',4]=_se[rwi]
		
		local row=`row'+1
	}	
	


	
clear
			
		svmat2 countryfigure

		drop if countryfigure1==.
		
		g country_code=_n
			run "$input/labels.do"
		
		g HI95_mean=(1.96*countryfigure2)+countryfigure1
		g LI95_mean=countryfigure1-(1.96*countryfigure2)

		g HI99_mean=(2.06*countryfigure2)+countryfigure1
		g LI99_mean=countryfigure1-(2.06*countryfigure2)
		
		
		g HI95_mean_ss=(1.96*countryfigure4)+countryfigure3
		g LI95_mean_ss=countryfigure3-(1.96*countryfigure4)

		g HI99_mean_ss=(2.06*countryfigure4)+countryfigure3
		g LI99_mean_ss=countryfigure3-(2.06*countryfigure4)
		
				
		set scheme plotplain
		
		
		rename (countryfigure1 countryfigure3 countryfigure5 countryfigure7 countryfigure9 ) (mean_beta mean_ss_beta max_beta urban_beta rural_beta)
		
		sort mean_beta
		g n=_n
		g n2=n+0.1
		
		

		label val country_code country_code
			decode country_code, g(country)
			replace country="IND/PAK" if country=="IND_PAK"
			labmask n, values(country)
			
		
		
		export delimited "$path/figdat/country_coefficients_collapsed.csv", replace
		
		twoway (rspike HI99_mean_ss LI99_mean_ss n if n<36, lc(gs10) horizontal) ///
					(scatter n mean_ss_beta  if mean_ss_beta<-3.9 & n<36, mc(cranberry%90) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>-4 & mean_ss_beta<-2.9& n<36, mc(cranberry%70) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>-3 &  mean_ss_beta<-1.9& n<36, mc(cranberry%50) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>-2  & mean_ss_beta<-0.9& n<36, mc(cranberry%30) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>-1 &  mean_ss_beta<0& n<36, mc(cranberry%10) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>0  & mean_ss_beta<1& n<36, mc(navy%10) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>1  & mean_ss_beta<2& n<36, mc(navy%30) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>2  & mean_ss_beta<3& n<36, mc(navy%50) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>3  & mean_ss_beta<4& n<36, mc(navy%70) ms(S) msize(large)) ///
					(scatter n mean_ss_beta  if mean_ss_beta>4& n<36, mc(navy%90) ms(S) msize(large) yti("") ylabel(, alternate valuelabel angle(0) labsize(large) nogrid) ysc(noline)) ///
					(rspike HI99_mean LI99_mean n if n<36, lc(gs10) horizontal) ///
					(scatter n mean_beta  if mean_beta<-3.9 & n<36, mc(cranberry%90) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-4 & mean_beta<-2.9& n<36, mc(cranberry%70) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-3 &  mean_beta<-1.9& n<36, mc(cranberry%50) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-2  & mean_beta<-0.9& n<36, mc(cranberry%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-1 &  mean_beta<0& n<36, mc(cranberry%10) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>0  & mean_beta<1& n<36, mc(navy%10) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>1  & mean_beta<2& n<36, mc(navy%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>2  & mean_beta<3& n<36, mc(navy%50) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>3  & mean_beta<4& n<36, mc(navy%70) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>4& n<36, mc(navy%90) ms(D) msize(large) ///
					yti("") ysc(noline axis(1)) ylabel(,nolab alternate valuelabel angle(0) labsize(large) nogrid noticks) ysc(noline) xlabel(,nogrid labsize(large))  xtitle("Country level correlation between" "wealth and average pollution",size(vlarge)) ysize(2.5) xsize(5) legend(off) xline(0))
					
					graph export "$figures/country_figure_mean_collapsed_ss_comp.png", as(png) replace
					
		
		twoway 	(rspike HI99_mean LI99_mean n , lc(gs10) horizontal) ///
					(scatter n mean_beta  if mean_beta<-3.9, mc(cranberry%90) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-4 & mean_beta<-2.9, mc(cranberry%70) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-3 &  mean_beta<-1.9, mc(cranberry%50) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-2  & mean_beta<-0.9, mc(cranberry%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>-1 &  mean_beta<0, mc(cranberry%10) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>0  & mean_beta<1, mc(navy%10) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>1  & mean_beta<2, mc(navy%30) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>2  & mean_beta<3, mc(navy%50) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>3  & mean_beta<4, mc(navy%70) ms(D) msize(large)) ///
					(scatter n mean_beta  if mean_beta>4, mc(navy%90) ms(D) msize(large) ///
					yti("") ysc(noline axis(1)) ylabel(,nolab alternate valuelabel angle(0) labsize(large) nogrid noticks) ysc(noline)  xlabel(,nogrid labsize(large)) xtitle("Country level correlation between" "wealth and average pollution",size(vlarge)) ysize(13) xsize(3) legend(off) xline(0))
					
					
					graph export "$figures/country_figure_mean_collapsed.png", as(png) replace

	
clear
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	